{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial4: Convolutional Layers - Spectral methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Why convolution in ML\n",
    "- Some theory on convolution\n",
    "- Convolution on graphs\n",
    "- Spectral-convolutional layers in PyTorch Geometric"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams.update({'font.size': 16})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why convolution in ML"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Weight sharing\n",
    "- Detection of translational invariant and local features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](fig/fully_connected.png)\n",
    "[Source](https://missinglink.ai/guides/convolutional-neural-networks/fully-connected-layers-convolutional-neural-networks-complete-guide/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](fig/cnn.gif)\n",
    "[Source](https://commons.wikimedia.org/wiki/File:Convolutional_Neural_Network.gif)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](fig/Convolution.png)\n",
    "[Source](https://commons.wikimedia.org/wiki/File:Convolution.PNG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some theory on convolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](fig/Convolution_of_box_signal_with_itself2.gif)\n",
    "[Source](https://en.wikipedia.org/wiki/File:Convolution_of_box_signal_with_itself2.gif)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "c[n] = (v * w)[n] = \\sum_{m=0}^{N-1} v[m] \\cdot w[n-m]\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def conv(v, w):\n",
    "    c = np.zeros(v.shape)\n",
    "    for n in range(len(v)):\n",
    "        c[n] = 0\n",
    "        for m in range(len(v)):\n",
    "            c[n] += v[m] * w[n - m]  \n",
    "    return c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD9CAYAAACcJ53WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXRU15nv/e8uVWlCICGEJNBQGGNmzCQzGyRRwsZJcOcmcW7nxsm6uY7TeZO8ndtvkk466dvdt293kk5y01lZ6awYZ61O287gpJOYtBvsKk3MZgYzm6k0oAEQEpqlUu33j6IIxghVHVXVqVP1fNZiCUradX4cikdHT+29j9JaI4QQwjpsZgcQQggRHincQghhMVK4hRDCYqRwCyGExUjhFkIIi7FH+wB5eXl6xowZhsf39vYyYcKEyAVKMnL+xkfO3/jI+TPu8OHD17XWU+/3uagX7hkzZnDo0CHD4+vq6igvL49coCQj52985PyNj5w/45RS3tE+J60SIYSwGCncQghhMVK4hRDCYgwVbqXUDqWUVkr9n0gHEkII8WBhF26l1J8Ci6OQRQghRAjCKtxKqRzg+8BfRCeOEEKIsYQ7HfCfgFNa618opX4ejUBCiIBj7cc41HaIsoIyluQvMTuOYYODg3R0dNDd3c3IyIjZcUyXmppKXl4e2dnZhp8j5MKtlFoHfAJpkwgRdcfaj/GpNz7FiH+E1JRUtm7aasnirZSioaGByZMnM2PGDBwOB0ops2OZRmtNf38/TU1NpKWlkZ6ebuh5QircSikH8BPgu1rrcyF8/fPA8wAFBQXU1dUZCgfQ09MzrvHJTs7f+Jh1/l67+RrD/mEAhkaGeHXfq3Rmd8Y8x3jZbDbS09NJS0tjcHCQwcFBsyPFhQkTJnDixAn6+voMjQ/1ivsvgQzgH0L5Yq31C8ALAGVlZXo8K6dk5dX4yPkbH7PO34EDB+BW4PeOFAfPrH7GklfcR44coaCggNTUVLOjxJX09HT6+vpYsWKFofFjFm6lVCnwdeA5IE0plXbXp9Nuv2HZrbWW5pUQEXKm4wx2ZcenfXx+6ectWbQh0CpxOBxmx4g7drsdn89neHwos0pmAunAy8DNu34BfOn27xcZTiCEeJcb/Tc40n6ETyz4BFmOLC51XjI70rgkc097NOM9J6G0So4BFfd5vJZAMf8pcGFcKYQQd9Q01uDXfp566Cna+tqobazF5/dht0V9TzhhEWO+ErTWnUDdvY/f/o7h1Vq/53NCCOM8Xg8lE0uYPXk2VaVVvH7pdQ61HWLVtFVmRxNxQvYqESKOdA12caDlAC6nC6UUa4rWkGHPwOP1mB1NxBHDhVtrrbTW34hkGCGSXV1jHT7to6q0CoAMewbritZR3VCNX/tNTifihVxxCxFHPF4PhRMKWZi38M5jVc4qrvdf51j7MROTiXgihVuIONE73Mveq3txlbreNetgffF6Um2puL1uE9PFl8Pem/yo9gKHvTfH/uIIevXVV1FKceLEifd8bvPmzSxZEptpm/I2tRBxYmfTTob8Q7icrnc9PsExgTXT11DdUM1XHvtKQkyv+7s/nOL01VuGxnYPDHO2tRu/BpuCuYUTmZge/lzx+dMn8TcfWBDWmC1btpCdnc3LL7/MP/3TP915vK2tDY/Hw7e+9a2wcxghV9xCxAm3182U9CksmfreqzaX00VLbwunbpwyIVl8uTXgw68Dv/frwJ9jJT09nY985CP8/Oc/x+//43sOv/jFL9Ba87GPfSwmOeSKW4g40O/rZ3fzbj4w8wOk2FLe8/nyknLsyo7b635X/9uqwr3Svdth703+24v7Gfb5cdht/OC/LmW5c3IE0z3Ys88+y4svvkhNTQ0uV+Cno5deegmXy8W0adNikkGuuIWIA3ub99Lv639PmyQoOy2bFdNW4PF60FrHOF18We6czCvPreIvNs3hledWxbRoAzz++OPMmDGDl156CYAzZ85w5MgRnn322ZhlkMItRBxwN7jJTsumrLBs1K9xOV00dDdw/ub5GCaLT8udk/lcxayYF20ILD78+Mc/zm9/+1v6+vp46aWXyMrK4oMf/GDMMkjhFsJkQyND1DfWU1FSgcM2+ptslSWV2JQNT4MsxjHbs88+S09PD7/97W955ZVX+NCHPkRmZmbMji+FWwiT7W/ZT89wD1XOqgd+3ZSMKSzLXyarKOPA7NmzWblyJV/96ldpaGiIaZsEpHALYTqP10OWIyukvUhcThcXOi9wuetyDJKJB3n22Wdpbm6mqKiIior77cMXPVK4hTCRz++jtrE2sMgmZeybDbhKA29eylW3+T73uc+htaapqQmbLbalVAq3ECY61HaIzsHOMdskQQUTCnh06qOyijLJSeEWwkQer4cMewZri9aGPKaqtIozHWdo6m6KYjIRz6RwC2ESv/ZT3VDNuqJ1ZNgzQh4XnOtd3VAdrWgizknhFsIkx9qPcb3/+p2+daiKJxYzL3eetEuSmBRuIUzi9rpx2BysL14f9liX08Xxa8dp622LQjIR76RwC2ECrTXVDdWsmb6GrNSssMdLuyS5SeEWwgSnbpyipbdl1L1JxjIzeyYPZz8sqyiTlBRuIUzg9rqxKzsVJcYXbricLg63HaZjoCOCyYQVSOEWIsa01ni8Hh4rfIzstGzDz1PlrMKv/dQ01EQwnbACKdxCxNj5m+dp6G4w3CYJmj15NiUTS2QVZRKSwi1EjHkaPCgUlaWV43oepRQup4u3Wt6ia7ArQumEFUjhFiLGPF4PywqWkZeRN+7nqiqtwqd91DfVRyCZsAop3ELE0OWuy1zovBDy3iRjWZi3kMIJhbIYJ8lI4RYihoL96I2lGyPyfEopXKUu9jbvpXe4NyLPaQmNB2DX9wIfY+jQoUMopdi9e/edx374wx+ilOIb3/jGncfeeecdlFL853/+Z1RyyM2ChYght9fNo3mPUjihMGLP6XK6ePnMy+xq2sWTDz0ZseeNqu1fhda3jY0dvAVtJ0H7QdmgYCGkTQr/eQoXweZvhTVk2bJl5OTkUFNTw7p16wCoqakhIyODmpo/zu6pqakhJSWFxx9/PPxcIZArbiFipKm7iTMdZ8Y9m+ReS6YuYUr6lORplwx0BYo2BD4OxO6NWZvNxvr166mtrQXA7/dTX1/PZz/7WQ4ePEhPTw8AtbW1lJWVMXHixKjkkCtuIWIkuDw90oU7xZbCxtKN/OHSHxjwDZBuT4/o80dFmFe679J4AH62BUaGICUVPvQilKyIXLYxVFRU8NWvfpWBgQFOnz5NZ2cnX/nKV/jJT37Crl272Lx5M3V1dXzqU5+KWga54hYiRtxeN3Nz51IysSTiz+1yuuj39bPn6p6IP3fcKVkBn9wGlV8PfIxh0QaorKxkcHCQvXv3Ultby+LFiykoKGDdunXU1tZy6tQp2traono7M7niFiIG2nrbOH7tOJ9f8vmoPH9ZYRnZadl4vJ6IvfEZ10pWxLxgBy1atIi8vDxqamo4evQolZWB+fiVlZW8+uqrlJSUkJqaytq1od8cI1xyxS1EDATbJJGaBngvh81BRUkF9Y31DI8MR+UYIkApxYYNG3C73ezatetdhfvo0aP87ne/Y+XKlWRmZkYtgxRuIWLA0+BhZvZMZubMjNoxqpxVdA93s79lf9SOIQIqKys5cOAAfX19d2aOLFu2jEmTJlFbWxv1u75L4RYiyjoGOjjcdjjqLYxV01YxwTFBtnqNgWBhLisrY9KkwFTE4IyTuz8fLdLjFiLKahpq8Gt/1NokQakpqawvXk9NQw1/veqvsdvkv3e0zJs3D631ex5/7bXXYnJ8ueIWIso8Xg9FWUXMzZ0b9WNVOavoHOzkcNvhqB9LmEcKtxBR1DXYxVstb1HlrEIpFfXjrZ2+lvSU9ORZjJOkpHALEUX1TfX4tC/ii25Gk+nIZF3RujvtGZGYxizcSqknlFI1SqlWpdSgUqpJKfWqUmp+LAIKYWVur5v8zHwW5S2K2TFdThfX+q9x/NrxmB1TxFYoV9y5wGHg88Am4GvAAmC/UsoZxWxCWFrvcC97m/fiKnVhU7H74XZD8QYcNoe0SxLYmK8mrfUvtNZf1lr/Rmtdr7V+CfgvwETgw1FPKIRF7WraxZB/KGZtkqCs1CxWT19Ntbf6vjMfhPUZvQy4cfujLNESYhS/Pv9rMuwZpKiUmB/bVeriau9V/vGtf+RY+7GYH19EV8iFWymVopRKVUo9AvwEaAV+GbVkQljYgZYDHGg9QL+vn8+4PxPz4pmfmQ/Ar879ik+/+Wkp3gkmnBn6bwHLb//+AlCptW6/3xcqpZ4HngcoKCigrq7OcMCenp5xjU92cv7Gx+j5+9n1n935/dDIEK/ue5XO7M4IJnuwN7veBECjTTl+0KRJk+ju7o75ca1gYGDA8P/NcAr3s8AkYCbwJcCtlFqntb5y7xdqrV8AXgAoKyvT5eXlhsIB1NXVMZ7xyU7O3/gYPX+/dP8SeiFFpeCwOXhm9TMsyV8S+YCjyGnPYfsb2/H5fdht9pgfP+jo0aNRu5mA1aWnp7N06VJDY0NulWitz2it39Ja/wLYCGQBXzV0VCES2PDIMCeunWBd0To+v/TzbN20NeZFc0n+Ev65/J8BeN/M95lStEX0GNrMQGvdqZS6AMyKcB4hLG9/y366h7v56JyPUl5SblqODSUbWJa/jLevG7y3o4hbhmaVKKUKgLnAxcjGEcL6PA0eMu2ZrJ6+2uwoVDmruNB5gStdV8yOIiIolJWTv1NK/bVS6mmlVIVS6jNAPeADvhf1hEJYiM/vo6ahhg3FG0hLSTM7zp055Im21eux9mO8+PaLps2WOX78OB/84AeZMmUKGRkZzJkzh29+85sxO34orZL9wDPA/wekAo1AHfDN+70xKUQyO9x2mM7BzpgvuhlN4YRCFuUtwu1189yi58yOc8e3D3ybsx1nDY3tGerh3M1zaDQKxZzJc8hKzQr7eebmzuUvV/xl2OMOHDhAeXk5s2bN4vvf/z7FxcW88847nDhxIuznMmrMwq21/jbw7RhkEcLy3F436SnprCtaZ3aUO1xOF98//H2ae5opyioyO864dQ93owmsCNVouoe7DRVuo770pS8xZcoU9u/ff+f2ZMHbl8WK7LQuRIT4tZ+ahhrWFq0l0xG9+w2Gq6q0iu8f/j4er4dPLvik2XEADF3pBh1rP8an3/w0w/5hHDYH33r8WzGbNdPX18eePXv48pe/HNV7So5FCrcQEXL82nGu9V+LmzZJUMmkEuZMnhNXhXs8luQvYeumrRxqO0RZQVlMpzrevHkTv99PcXFxzI55P1K4hYgQt9eN3WZnQ/EGs6O8h8vp4kfHfkR7X/ud5fBWtiR/iSlz0ydPnozNZqO5uTnmx76b3EhBiAjQWlPtrWb1tNVMTI2/lYLB+11WN1SbnMTaMjMzWbduHS+//DL9/f2m5ZDCLUQEnL5xmqu9V6N+Q2CjHs55mIeyH8LjTaxpgWb47ne/y40bN1i9ejUvvfQStbW1/PSnP+ULX/hCzDJI4RYiAtxeNykqhYqSCrOjjMpV6uJQ2yE6BjrMjmJpjz32GHv27KGkpIQvfOELPPXUU3znO9+Jad9betxCjJPWGk+Dh7LCMnLSc8yOM6oqZxVb395KbUMtH5r9IbPjWNrSpUv5wx/+YNrx5YpbiHF6p/MdvLe8VJXGZ5skaG7uXIqyinA3yC3NrE4KtxDj5PF6UCg2OjeaHeWBlFJUOat4q+Utbg3dMjuOGAcp3EKMk9vrZmn+UvIy8syOMiaX04XP76O+sd7sKGIcpHALMQ5Xuq5wofNC3C26Gc2ivEXkZ+bLHeAtTgq3EOMQ3HXPVWqNwm1TNlylLvZe3UvfcF9Mjil3mn+v8Z4TKdxCjIPb62bhlIVMy5pmdpSQuZwuBkcG2dm8M+rH8vv9pi5UiVf9/f04HA7D46VwC2FQc08zp2+ctkybJGhZ/jJy03Njshinr6+PpqYmOjo6GB4eTvqrb601fX19NDc3k59vfOsBmccthEHBwhevqyVHk2JLobK0ktcvvc6Ab4B0e3rUjjUyMkJpaSnXrl3jxo0b+Hy+qB3LKhwOBwUFBUyaNMnwc0jhFsKg6oZqZk+eTemkUrOjhK2qtIrfnP8N+67uo6I0uqs909PTKSkpieoxko20SoQw4FrfNY61H7NcmyTosWmPMSl1UsLd0ixZSOEWwoDqhmo0Ou5XS47GYXNQXlJObWMtwyPDZscRYZLCLYQBHq+HGZNm8HDOw2ZHMazKWUX3UDcHWg+YHUWESQq3EGG6OXCTQ22HqHJWoZQyO45hq6evJtOeKYtxLEgKtxBhqm2sZUSPWLa/HZSWksaG4g2Bv49/xOw4IgxSuIUIk9vrpiiriHm588yOMm4up4uOgQ6OtB8xO4oIgxRuIcJwa+gW+1v24yp1WbpNErSuaB3pKenSLrEYKdxChKG+sR6f32f5NklQpiOTtUVrqfZW49d+s+OIEEnhFiIMHq+H/Ix8Hp36qNlRIsbldNHe386JayfMjiJCJIVbiBD1Dfex5+oeNjo3YlOJ819nQ/EG7Da73EjYQhLn1SdElO1q3sXgyKDl9iYZy8TUiayethpPgyfpN4GyCincQoTI4/WQm57LsvxlZkeJuCpnFc09zZzpOGN2FBECKdxChGBwZJCdTTupKKkgxZZidpyIqyipIEWlSLvEIqRwCxGCvc176fP1JVybJCgnPYeywjLcXre0SyxACrcQIfA0eJiYOpEVhSvMjhI1VaVVXLl1hYudF82OIsYghVuIMQyPDFPbWEtFSQWOFOO3m4p3G50bUSjcDbIYJ95J4RZiDAdaD9A91G2ZGwIblZeRx9L8pdLntgAp3EKMwe11k2nPZE3RGrOjRJ3L6eL8zfM03GowO4p4ACncQjyAX/upbaxlffF60lLSzI4TdcGfKmTvkvgmhVuIB7g4eJGOgY6E2ZtkLNOyprFwykJpl8S5MQu3UurDSql/V0p5lVL9SqlzSqlvKqUmxiKgEGY61neMtJQ0Hi963OwoMeNyujh54yQtPS1mRxGjCOWK+0vACPBXwJPAj4HPAm6lEmjDBiHu4dd+jvcdZ+30tWQ6Ms2OEzPBuepyI+H4ZQ/haz6gtb5215/rlVIdwM+AcqAmGsGEMNuJayfoGulKmjZJUOmkUmZPno3H6+HZ+c+aHUfcx5hXzPcU7aCDtz8WRTaOEPHj52d+jkIxJWOK2VFizuV0caT9CD848gOOtR8zO464h9FWx4bbH2VHGpGQjrUfY8eVHWg0f17z50lXvEqySgD46ds/5dNvfjrp/v7xLpRWybsopYqA/w14tNaHRvma54HnAQoKCqirqzMcsKenZ1zjk52cP2N+fePXaAJ7dgyNDPHqvlfpzO40OVXs7OncA4BGj+vvL6+/6AircCulsoDXAB/w30f7Oq31C8ALAGVlZbq8vNxwwLq6OsYzPtnJ+TPGvcsNPaBQpKak8szqZ1iSv8TsWDGT057D9h3bGdEjOFIchv/+8vqLjpBbJUqpdGAbMBN4QmvdFLVUQphIa82J6yeYnzuf9+e8n62btiZV0QZYkr+Ev13ztwA8O//ZpPv7x7uQCrdSygH8O7ACeEpr/XZUUwlhooudF7ly6woffOSDbMrelLRF6+mHn6Yoq4hzHefMjiLuEcoCHBvwCrAReFprvT/qqYQwkbvBjUKxsXSj2VFMpZTCVepiX8s+uoe6zY4j7hLKFfePgI8A3wV6lVKr7vpVHN14QsSex+th8dTFTM2canYU07mcLnx+H/VN9WZHEXcJpXBvvv3x68C+e349F6VcQpii4VYD52+eT7pFN6N5dOqj5Gfky94lcWbMWSVa6xkxyCFEXAjuiieFO8CmbFSWVvL7C7+nb7gvqZb+xzPZa0SIu3i8HuZPmU9RliwKDqpyVjEwMsDu5t1mRxG3SeEW4raWnhZO3jiZsDcENmpZwTImp02WdkkckcItxG3B3fAS/RZl4bLb7FSWVlLfVM/gyKDZcQRSuIW4w+P1MCtnFjOyZ5gdJe64nC76fH3su7rP7CgCKdxCAHC9/zpH249Km2QUKwtXMtExUW5pFiekcAsBVHur0WiZTTIKR4qD8pJy6hrrGPYPmx0n6UnhFoLAaknnJCeP5DxidpS45XK6uDV0i4MtB8f+YhFVUrhF0usc6ORQ6yFcpS6UUmbHiVtrpq8hw56Bu0HaJWaTwi2SXm1jLSN6RPrbY0i3p7O+eD01DTWM+EfMjpPUpHCLpOf2upk+YTrzp8w3O0rcczlddAx0cKT9iNlRkpoUbpHUuoe62d+yn43OjdImCcH6ovWkpaTJYhyTSeEWSW1n006G/cPSJglRpiOTNdPX4Gnw4Nd+s+MkLSncIql5vB6mZkxl8dTFZkexjCpnFe197bx9Xe6nYhYp3CJp9Q33sbt5N5WlldiU/FcI1YaSDdhtdmmXmEherSJp7bm6h4GRAWmThGlS6iRWTluJ2+tGa212nKQkhVskLbfXTU5aDssLlpsdxXKqSqto7mnmbMdZs6MkJSncIikNjQyxs2knlaWV2G1j3k9E3KOitAKbssneJSaRwi2S0r6r++gd7pUtXA3KTc+lrKDszla4IrakcIuk5Pa6meiYyKppq8yOYlkup4vLXZe52HnR7ChJRwq3SDrD/mFqG2vZULIBR4rD7DiWtbF0I4C0S0wghVsknYOtB7k1dEu2cB2n/Mx8lkxdItMCTSCFWyQdj9dDhj2DtdPXmh3F8lxOF+dunqPxVqPZUZKKFG6RVEb8I1Q3VPN40eOk29PNjmN5wZ9aZKvX2JLCLZLK0fajdAx0yKKbCCnKKmL+lPnSLokxKdwiqXgaPKTaUnm8+HGzoySMKmcVb19/m9beVrOjJA0p3CJp+LUfj9fDmqI1THBMMDtOwgjOhZer7tiRwi2SxsnrJ2nra5M2SYTNyJ7BrJxZMi0whqRwi6Th8XqwKzsbijeYHSXhVDmrONp+lOv9182OkhSkcIukoLXG7XWzctpKstOyzY6TcFxOFxpNTUON2VGSghRukRTO3TxHU0+TLLqJkkdyHsE5ySntkhiRwi2SgtvrxqZsVJZWmh0lISmlcJW6ONh6kM6BTrPjJDwp3CIpeLwelhcsJzc91+woCavKWcWIHqG2sdbsKAlPCrdIeJc6L3Gp65Js4Rpl86fMZ/qE6bLVawxI4RYJL9h3De5mJ6JDKcVG50b2Xd1Hz1CP2XESmhRukfA8DR4WT11MwYQCs6MkvCpnFcP+Yeqb6s2OktCkcIuE1tjdyNmOs7LoJkYWT13M1IypsooyykIq3EqpYqXUD5VS+5RSfUoprZSaEd1oQoxfsIBImyQ2gjN3djfvpm+4z+w4CSvUK+5ZwDPATWBX9OIIEVker4d5ufMonlhsdpSkUeWsYmBkgL1X95odJWGFWrh3aq0LtNZPAb+OZiAhIqW1t5UT109ImyTGlhcsJyctRxbjRFFIhVtr7Y92kKhoPAC7vhf4mIzjLe6w9yY/qr3AYe9NQ+O/s28rAEMDkw0f/z8uDhk+/njzW5XdZqeytJLahlq2d27nWPsxsyMlHLvZAaKm8QD86/thZBhS7LDpHyB/Xujj28/Am1+HEZ+J478Bfh+kpMInt0HJitDHW9xh700+tnU/wyN+7DYbX3/fPGYXTAx5fN2Vg7zRFPjh8Men/pFb3VmUz3gs5PHn27r5h9fPMDziZ9ul/WEfPzje5/eTarfxynOrWO409g3EimZmz6R/pJ/tXdupebOGrZu2siR/idmxEobSWoc3QKnngK3AQ1rrK6N8zfPA8wAFBQXLf/nLXxoO2NPTQ1ZWVtjjHjn3Y4padhg+bjzxY+PKQ/+NBueHwx5r9PyZ7dfnBnn9ss/w+NS8N0jNq0Up0FoxdG0TQzcqIpgwdDbgvzzi4P0Pp5pyfDPs6NzB612vA2DDxvty3sem7E0mp7KWioqKw1rrsvt9LipX3FrrF4AXAMrKynR5ebnh56qrq8PQ+Gv/Bi2AsoFtHFfMfp9547d/BdDY7GnMrPwEMw1ccRs+fybb3XsaLl/GpjB0xb311EEO3gwUbbSd/77MZfiK25ES/vHPt3Xzd384hV9DqsPGn7oeS6or7pz2HLbv2I5f+3GkOHhm9TNyxR1Bidkq0RqaDkLxSpjzBMx4PPw2w0OPw/QlcGWXeeObD8OJX8FHX06qNgnAqeZbFE3O4GMrSlk1c0rYRe9fL7UwuTePRzKeYNPMtXz00fBuVbb64SksLMrmF56Dhoru6oen0NrVz4/rL/H3Ty9MqqINsCR/CV9Y+gV+cOQHfHHZF6VoR1hiFu6rR6GrAcr/EpZ+3PjzlKwYX8Ec7/hVfwYnfgk9yXUvvxs9g7x1+Qafq5jF5ypmhT2+a7CLAy0H+MSCT/A/l/9PwzmWOyfT/XCq4aL73OMz+cnOS3hvJOd85o/P+zg/PvpjLnddNjtKwknMlZNntoFKgTlPmZ1kfKYtgZxSOL3N7CQx5T7dhl/DkwsLDY2va6zDp32mTwOckpXGyoemsP1ki6k5zJJuT2dBxgKqG6oZ8Y+YHSehhFy4lVIfVkp9GFh++6HNtx+Lr/tAaR0odA89DpkW38JTKZi3BS7VwsAts9PEzPaTrZTmZjJ/2iRD4z1eD9MmTGPBlAURTha+zYsKuXitl3faus2OYorFmYu5MXCDY9dkSmAkhXPF/evbv/7s9p//5faf/y7Socal/TR0XAwUvEQwbwuMDMH5N8xOEhNd/cPsvXidzQsLUUqFPb53uJe9V/eysXSjofGR9sSCwE8N208mV7sraEHGAlJtqbJ3SYSFXLi11mqUX+VRzBe+09sABXPfb3aSyCh+DCZOgzOvmZ0kJqrPtDE8og23SXY27WTIP2R6mySoYFI6y52Tk7Zwp9vSWVO0Bk+Dh3CnHovRJV6P+8w2KF0NExNkC0+bLfBN6B0PDPWanSbqtp9sZVp2OouLcwyNd3vd5GXkxdUshs0LCznTcgvvjcT/97ufKmcVrb2tnLx+0uwoCSOxCvf1C4FWyfwEaZMEzd8Cvn64kNg/bvYO+lDNGIAAABMRSURBVNh5/hpPLCjEZgu/zdHv62d38242lm7EpuLnpZ3s7ZINxRuwKzvuBtm7JFLi59UdCcF2wrwPmJsj0krXQOaUhJ9dUnuunUGfn80G2yR7m/fS7+uPuzu5l+RmsqgoO2kLd3ZaNiunrcTjlXZJpCRW4T69DYqWQ3aCbeGZYoe57wu8QekbNDtN1Gw/2UpeViplM4zNBnI3uMlJy6Gs4L6rhE315MJCjjd2crWz3+wopnA5XTR2N3L+5nmzoySExCncN73QcixxZpPca97TMNQNFxPzDtoDwyPUnm1n04JCUgy0SYZGhqhvrKeipAK7Lf7WlQV/itiRpFfdlaWV2JRNtnqNkMQp3Gf+EPiYaP3toIfWQ1p24M3XBLTz/DX6hkYMt0n2t+ynZ7gn7tokQTOnZjGnYGLSFu7c9FyWFyyXaYERkkCFexsULILcmWYniQ57KszZDGdfD2xVm2B2nGwlO8PBqplTDI33eD1kObJYNW1VhJNFzpMLCzno7aC9e8DsKKZwlbq42HWRS12XzI5ieYlRuG+1QONbiXu1HTR/Cwx0BjauSiBDPj/uM21UzS/AkRL+S9Ln91HbWMuGkg2kpsTv1qmbFxWiNbx5qs3sKKYI3vdTrrrHLzEK99n/CHxM1P520MOV4JiQcLNL9l68TveAz3Cb5FDbIToHO6kqjY9FN6OZUzCRh/ImJG27pGBCAYunLpbCHQGJUbhPvwZ5syF/rtlJosuRAbM3Bb5RJdCmPTtOtpKVZmfdI3mGxnu8HjLsGawpWhPhZJGllOLJhYXsu3SDm71DZscxRZWzijMdZ2jsbjQ7iqVZv3D3XgfvnsS/2g6atwV6r0HDfrOTRIRvxM+bp9uonJtPmj0l7PF+7ae6oZp1RevIsGdEIWFkbV5YyIhf4z6T3O2Sam+1yUmszfqF++zroP2J398OemQT2NMTZnbJgSsddPQOGW6THGs/xvX+63GzN8lYFhVlU5STkbTtkuKJxczLnSerKMfJ+oX7zDbIcULho2YniY20LHh4Y2D6o99vdppx23GylXSHjQ1zphoa7/a6SbWlsr54fYSTRUewXbL7net0DyTe7KBQVDmrOHHtBK29yfnNKxKsXbj7O+FSfeBqOw628IyZ+VvgVjNcPWJ2knHx+zU7TrZSPjufzNTwF81oraluqGbN9DVMcEyIQsLo2LywkKERPzVn282OYorgXPvqBmmXGGXtwn1+B/iHA6sKk8nsJwM3ID5t7a1ejzbepL17kM2LjLVJTt04RUtvS9wuuhnNstLJTJ2Yxva3k/OK86Hsh5iVM0tml4yDtQv36W0wcXpgf5JkkpEDD20ItIksvGnP9rdbcaQoKubmGxrv9rqxKzvlJeWRDRZlNpviiQUF1J1vp2/IZ3YcU7icLo60H+FG/w2zo1iSdQv3YA9crA7sBGiz7l/DsPlb4OYVaH3b7CSGaK3ZfrKVdbPymJTuMDTe4/WwYtoKstOyo5AwujYvnMbAsJ+d56+ZHcUUrlIXfu2nprHG7CiWZN2K986b4BtIntkk95r7flA2y84uOdl8i+bOfjYvnGZo/Pmb52nobrBcmyRo5UO5TM50JO1Wr7Mnz6Z0Yqm0SwyybuE+sw0mTA3c7SYZTcgD51rLrqLcfrKFFJuiar6xOxV5GjzYlI3KksoIJ4sNe4qNqvkF1JxpZ9CXOIupQqWUwuV0caDlAF2DXWbHsRxrFu7hfjj/ZmCPalv4izYSxrwtcP0cXDtndpKwaB2YTbJqZi6TJxjbW8Tj9bAsfxlTMoxtShUPNi+cRvegjz0XrpsdxRRVzip82kddY53ZUSzHmoX7Yg0M9ybPasnRzLt9Q2SLXXWfb+vh0vVenjTYJrncdZkLnRcs2yYJWjNrChPT7Ek7u2TBlAUUTiiUdokB1izcp7dBek5gj+pkNmk6FK+w3B3gt59sQSl4YoHBNsnt/+jB5dNWlWZPYeO8fNxn2hgesf5iqnAppXCVuth7dS+9w8l5I2WjrFe4fUNwbjvMeQpSwp+NkHDmbwnMLOm4bHaSkO042UqZczL5E9MNjXd73Tya9yiFE4zN/44nTy6cRmffMG9d6jA7iilcThdD/iF2Nu00O4qlWK9wX94Jg13JO5vkXsEbI1tkdsnl672cbe023CZp6m7iTMcZy7dJgjbMnkqGI4XtJ1vMjmKKJVOXMCV9itzSLEzWK9xnXoPULJhZYXaS+DB5BkxbbJk+d7BAPWlwU6ngMulEKdwZqSlUzJ3KG6faGPFbdzGVUSm2FDaWbmR38276fcl5I2UjrFW4R3yB3QBnPwEOYz9mJ6R5W6D5EHQ1m51kTDtOtrK4OLBDnhFur5u5uXMpmVgS4WTmeXLhNK73DHLYe9PsKKZwOV30+/rZ27zX7CiWYa3C3bAX+m7IbJJ7zb+9V0vwhslxqulmHyeaugy3Sdp62zh+7Tiu0sS42g6qnJtPqt2WtO2SssIystOyZavXMFircJ/eBvYMeMQaey/HTN4jMHVe3Pe5g3tQG917O9gmscre26HKSrOz/pE83jjZirbw3jNGOWwOKkoqqG+sZ2gkOe8MFC7rFG6/P3BFOWsjpFpnC8+Ymb8FvHuhJ363Ct1xspW5hROZkWfs38/T4GFm9kxm5syMcDLzPblwGle7BjjelJyrCKucVfQM97C/JTHu7BRt1incTQehp/WPbQHxbvO2APqPN06OM+23BjjccNPw3iQdAx0cbjucMG9K3qtqXgF2m0radsmqaavIcmTJYpwQWadwn9kGNkfgjUnxXgULIHdm3M4ueeNUK1pjeO/tmoYa/NqfcG2SoOxMB6sfnsKOJG2XpKYE7mJU21iLz5+cW92GwxqFW+tAQXq4AtKtt4VnTCgVuOq+sgv64m8xx/aTrcycOoFH8rMMjfd4PRRnFTNn8pwIJ4sfmxdOw3ujjzMt3WZHMUWVs4rOwU4OtR0yO0rcs0bhbjkGXQ0ym2Qs87eA3xdYWRpHOnqHeOtyB5sXFqIM3GKua7CLt1reospZZWi8VWxaUIBNwY4kbZesLVpLhj1D2iUhsEbhPr0NVEpgN0AxuunLILsk7maXuE+3MuLXhvvb9U31+LQvYfvbQXlZaTw2Izdp9+jOsGewrmgd1Q3V+HXy7d0Sjvgv3FoHCtGMdZCZa3aa+KZUYAn8xRoYuGV2mju2n2yleHIGC6ZPMjTe7XVTkFnAwryFEU4WfzYvLOSd9h4utPeYHcUUrlIX1/uvc6z9mNlR4lrcF+4JvQ1w44LsTRKqeVtgZChwh6A40NU/zJ4L1w23SXqHe9nbvBeX04VNxf3LddyCi5OStV2yvng9DptD9i4ZQ9z/T8i7vg9QMPcDZkexhpKVkFUQN3eArznbxvCINrxaclfTLob8Qwm3WnI0hdnpLC3NSdp2SVZqFmumr6G6oTopZ9eEKqTCrZQqUUr9RinVpZS6pZT6rVKqNNrhAKZe2wulq2Cisb2bk47NFrgf5QUPDPWZnYbtb7dSMCmNpSU5hsa7vW5y03NZmr80wsni1+aFhZy6eouGG+b/+5nB5XTR0tvCqRunzI4St8Ys3EqpTKAGmAt8EngWeASoVUpFdwnjjYtk9XplNkm45m+B4b5A8TZR76CP+vPXeHJBITZb+G2SAd8Au5p3sbF0IylJdIu64Ju4O04lZ7ukoqQCu7JLu+QBQrni/jQwE/gTrfXvtdavAVsAJ/CZaIa7MztinrRJwuJcBxm5ps8uqTt3jUGf33CbZM/VPfT7+hN+Nsm9SnIzWTB9UtK2S7LTsnms8DE8Xo+0S0YRSuHeAuzXWl8IPqC1vgzsAaK6/vyto//GDycX8NrxNwyN/9WJXfyP33+bX53YlVzjU+ww9ymOXHiD109/h9fqXzR0/NfqX+R//ewjhsdv3/1TyvP/hbbG3xka/+tjW0lTKaRev2hoPI0HYNf3Ah8Nji/1/mZc440ef/PCQlTjAepe/CpnDxr7yensQQ/7fvZXpo7vO/orQ+NdThcN3Q188d+eNu31a/b4B1FjfUdTSrUCr2mtP3PP4/8CfERrPfVB48vKyvShQ+GvhPq/v/kK/9rzn2hAAY6RSfgJvTMzwhD+lD+uILSN5JJC6HcUt/r4DDrpSem/c/7yfZpUHXq7Ykhp2u0qYuOnjUBGGLNCBtA02wKvzXSt2dpjY4ktM+TxDPdDZwMEE+SUgiOMPcBvj9do1DjGGz1+b283Gb1XUbcTtKqpDNnSQh6f6h+kUF+z7PijqUP8r8LA6z0eXr/jGZ+m4a9nfpGnNzwX8ngApdRhrXXZ/T5nD2F8LnC/Hd47gMmjHPB54HmAgoIC6urqQkt6l8vtB9CZgFJorcny+xgiP+TxfbRBSmBqs9ag/A7Sk2i8gy50CnfOX5pW5PlCLxzX7YGiH6nxqX4osoV+84sm3R/4eVAphoE9NphO6PP4M4ebmIBGARpN77Cmz2Gd8SNDvWSgb//7a/pJpd02PeTx+f5msPD482nNgXMXJ6/f8Yz3odl9+vdk61khjx9LKIUbApcN9xr124/W+gXgBQhccZeXl4cdrEt9nH2X/hkfGruGv5j9fFjfsX51Yhd/f/jP0fhA2/mrlX/DRx99PGnGv1b/In9/1/l7fnZ43/EjPf65MMcfO/lzPn3wHxlG49CwdsPXyF/4sZDH03gAfhaY065SUsn6rz8lq2RF2OP9vkFs9jTD440e/+xBDzn/8ac4tI9h7Ay/7wesfiz0Xv/Zgx4GLTy+vf5FUuPo9Tve8evm/wnlG8pDHj+WUFolbcDvY90qgcBffvfp37Nu/p+E/WMGBIrfm5f2smnmmrCKXqKMH+/5e63+RQ5feYPlM54wZfyxkz/n0KU3KJv5BEvCKdpBjQcCm27NeBzCKbp3jb9U82/MrPyE4fHjOf7Zgx5unq5h8vxK5oZR9OJp/KXd/87MdR8yNN7s15/Z4x/UKgmlcNcAqVrrdfc8Xnd7/IYHjR9P4Qaoq6vDyBW7CJDzNz5y/sZHzp9xDyrcobxbtA1YpZS6c9sRpdQMYO3tzwkhhIihUAr3VuAK8JpS6mml1BbgNaAR+EkUswkhhLiPMQu31roXqATOAy8BrwCXgUqtdXJuYSaEECYKaVaJ1roB+FCUswghhAhB3O8OKIQQ4t2kcAshhMWMOR1w3AdQ6hrgHcdT5AHXIxQnGcn5Gx85f+Mj588452jrZKJeuMdLKXVotLmMYmxy/sZHzt/4yPmLDmmVCCGExUjhFkIIi7FC4X7B7AAWJ+dvfOT8jY+cvyiI+x63EEKId7PCFbcQQoi7SOEWQgiLibvCrZQqUUr9RinVpZS6pZT6rVKq1OxcVqGUKldK6fv86jQ7W7xRShUrpX6olNqnlOq7fZ5m3OfrJiulXlRKXVdK9SqlPEqpRbFPHF9COX9KqRmjvB61UirHnOTWF+odcGJCKZUJ1ACDwCcJ3Hnn/wC1SqlHb294JULz/wIH7/qzz6wgcWwW8AxwGNgFbLr3C5RSisD2xQ8BXyBwG7+vEXhNLtFaN8UubtwZ8/zd5Zu8dxvo7ijlSnhxVbiBTwMzgTnBu8orpU4A7wCfAf6vidms5ozWer/ZIeLcTq11AYBS6jnuX3i2AOsI7IZZe/tr9xHYIfMrBL5BJqtQzl/QJXk9Rk68tUq2APuDRRtAa30Z2AM8bVoqkZC01v4QvmwLcDVYtG+P6wL+QJK/JkM8fyIK4q1wLwBO3ufxU8D8GGexuleUUiNKqRtKqZ/L+wSGPeg1WaqUyopxHqv6plLKd/u9q23yHsH4xFurJJdAD/FeHcDkGGexqi7ge0A9cAtYCvwVsE8ptVRr3W5mOAvKJXAHqHt13P44GZAbioxukMCdst4ErgFzCbwe9yqlVmitz5gZzqrirXBD4A3Je6mYp7AorfVR4OhdD9UrpXYCBwj0Y79hSjDrUshr0jCtdQvwZ3c9tEsptYPATyxfBz5uSjCLi7dWyU0CVzj3msz9r8RFCLTWRwjceu4xs7NYUAejvyZBXpdh01o3AruR16Nh8Va4TxHoKd5rPnA6xlkSzWhXjuLBHvSabJD7rhomr8dxiLfCvQ1YpZSaGXzg9oT+tbx3DqgIkVKqDJgNvGV2FgvaBhQppTYEH1BKTQI+gLwmDbn9Rvla5PVoWFxtMqWUmgAcB/oJ9GI18PfAROBRuboZm1LqFQJzjI8AnQTenPwa0Acs01rL3UjuopT68O3fbiTQi/1/CLyJdk1rXa+UshH4sb4E+DJ/XIDzKLD49o/9SSuE8/c9AheI+24/PofA+csGVmqtz8U+tfXFVeGGO9+Nvw9UEfhxqhr4otb6ipm5rEIp9TXgTwEnkAm0AtuBv7n9RpG4i1JqtP8A9Vrr8ttfkwt8F/gTIJ1AEfoLrfXxmISMY2OdP6XUp4DPElhlOZHAbcxqgL+Tom1c3BVuIYQQDxZvPW4hhBBjkMIthBAWI4VbCCEsRgq3EEJYjBRuIYSwGCncQghhMVK4hRDCYqRwCyGExfz/+rUNEeELfqQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 20\n",
    "v = np.zeros(N)\n",
    "v[8:12] = 1\n",
    "w = np.zeros(N)\n",
    "w[1:5] = 1\n",
    "c = conv(v, w)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.gca()\n",
    "ax.plot(v, '.-')\n",
    "ax.plot(w, '.-')\n",
    "ax.plot(c, '.-')\n",
    "ax.legend(['v', 'w', 'c'])\n",
    "ax.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fourier transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transformation $\\mathcal F: \\mathbb{R}^N \\to \\mathbb{R}^N$ with\n",
    "\n",
    "\\begin{align*}\n",
    "\\mathcal F^{-1}(\\mathcal F (v)) &= v\\\\\n",
    "\\mathcal F(v * w) &= \\mathcal F(v) \\cdot \\mathcal F(w).\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This implies\n",
    "\\begin{align*}\n",
    "v * w &= \\mathcal F^{-1}(\\mathcal F (v * w))\\\\\n",
    "&= \\mathcal F^{-1}(\\mathcal F(v) \\cdot \\mathcal F(w))\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4.09862852, 4.68214615, 4.70242855, 4.15529448, 3.82105494,\n",
       "       4.52767695, 4.31882443, 4.57647736, 4.95389019, 4.91595713,\n",
       "       4.35010702, 3.81069648, 4.29596588, 5.10180936, 4.27236417,\n",
       "       4.74873021, 4.94188403, 4.08582287, 5.03542755, 4.42217054])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v, w = np.random.rand(N), np.random.rand(N)\n",
    "conv(v, w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4.09862852, 4.68214615, 4.70242855, 4.15529448, 3.82105494,\n",
       "       4.52767695, 4.31882443, 4.57647736, 4.95389019, 4.91595713,\n",
       "       4.35010702, 3.81069648, 4.29596588, 5.10180936, 4.27236417,\n",
       "       4.74873021, 4.94188403, 4.08582287, 5.03542755, 4.42217054])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.fft import fft, ifft # Fast Fourier Transform / Inverse FFT\n",
    "np.abs(ifft(fft(v) * fft(w)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition of the Fourier transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Fourier transform can be computed as\n",
    "\n",
    "\\begin{align*}\n",
    "\\mathcal F(v) = U\\cdot v, \\;\\;\\mathcal F^{-1}(v) = \\frac{1}{N}\\ U^H \\cdot v\n",
    "\\end{align*}\n",
    "\n",
    "where the $N\\times N$ matrix $U$ is defined as\n",
    "\\begin{align*}\n",
    "\\\\\n",
    "U = \n",
    "\\begin{bmatrix}\n",
    "u_0(0) & u_1(0) & \\dots & u_{N-1}(0)\\\\\n",
    "u_0(1) & u_1(1) & \\dots & u_{N-1}(1)\\\\\n",
    "\\vdots & \\vdots& & \\vdots\\\\\n",
    "u_0(N-1) & u_1(N-1) & \\dots & u_{N-1}(N-1)\\\\\n",
    "\\end{bmatrix} \n",
    "\\end{align*}\n",
    "\n",
    "and $u_0, \\dots, u_{N-1}$ are functions defined as\n",
    "\n",
    "\\begin{align*}\n",
    "u_n(x)&:= \\cos\\left(2 \\pi \\frac{n}{N} x\\right) - i \\sin\\left(2 \\pi \\frac{n}{N} x\\right).\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matrix_U(N):\n",
    "    u = lambda n, N: np.cos(2 * np.pi / N * n * np.arange(N)) - 1j * np.sin(2 * np.pi / N * n * np.arange(N))\n",
    "    U = np.empty((N, 0))\n",
    "    for n in range(N):\n",
    "        U = np.c_[U, u(n, N)]\n",
    "    return U\n",
    "\n",
    "\n",
    "def fourier_transform(v):\n",
    "    N = v.shape[0]\n",
    "    U = matrix_U(N)\n",
    "    return U @ v\n",
    "\n",
    "\n",
    "def inverse_fourier_transform(v):\n",
    "    N = v.shape[0]\n",
    "    U = matrix_U(N)\n",
    "    return (U.conj().transpose() @ v) / N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.77635684e-15-0.00000000e+00j,  8.88178420e-16-2.22044605e-16j,\n",
       "        0.00000000e+00-4.44089210e-16j,  2.22044605e-16+1.11022302e-15j,\n",
       "        2.22044605e-16+7.77156117e-16j,  2.24993635e-15+1.11022302e-15j,\n",
       "       -1.33226763e-15-3.08086889e-15j, -1.99840144e-15+2.44249065e-15j,\n",
       "       -2.22044605e-15+8.88178420e-16j, -3.05311332e-15+1.27675648e-15j,\n",
       "        4.44089210e-16-6.32292965e-15j, -5.27355937e-16-8.88178420e-16j,\n",
       "       -8.88178420e-16-8.88178420e-16j, -3.77475828e-15-3.88578059e-15j,\n",
       "        4.44089210e-16-3.83026943e-15j,  7.57380270e-15+3.77475828e-15j,\n",
       "        2.22044605e-16+3.05311332e-15j, -4.44089210e-16-7.38298311e-15j,\n",
       "       -6.66133815e-16+7.99360578e-15j, -8.65973959e-15+8.88178420e-16j])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fft(v) - fourier_transform(v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.11022302e-16-0.00000000e+00j,  4.16333634e-17+6.93889390e-18j,\n",
       "        0.00000000e+00+2.08166817e-17j,  1.38777878e-17-5.55111512e-17j,\n",
       "        1.38777878e-17-3.98986399e-17j,  1.12431765e-16-5.55111512e-17j,\n",
       "       -6.93889390e-17+1.54390389e-16j, -9.71445147e-17-1.24900090e-16j,\n",
       "       -1.11022302e-16-4.16333634e-17j, -1.52655666e-16-6.24500451e-17j,\n",
       "        2.77555756e-17+3.16146482e-16j, -2.68882139e-17+4.51028104e-17j,\n",
       "       -4.16333634e-17+4.16333634e-17j, -1.94289029e-16+1.94289029e-16j,\n",
       "        2.77555756e-17+1.91686944e-16j,  3.78603399e-16-1.87350135e-16j,\n",
       "        1.38777878e-17-1.52655666e-16j, -2.08166817e-17+3.67761377e-16j,\n",
       "       -3.46944695e-17-3.98986399e-16j, -4.37150316e-16-4.16333634e-17j])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ifft(v) - inverse_fourier_transform(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Connection with the Laplacian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The functions $u_n$ (the columns of the Fourier transform matrix) are eigenvectors of the Laplacian:\n",
    "\n",
    "\\begin{align*}\n",
    "u_n(x)&:= \\cos\\left(2 \\pi \\frac{n}{N} x\\right) - i \\sin\\left(2 \\pi \\frac{n}{N} x\\right)\\\\\n",
    "\\Delta u_n(x)&:= \\left(-4 \\pi\\frac{n^2}{N^2}\\right) u_n(x)\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "v * w \n",
    "= U^H ((U  w) \\odot (U  v))\n",
    "\\end{align*}\n",
    "\n",
    "or if $g_w=\\mbox{diag}(U w)$ is  filter\n",
    "\\begin{align*}\n",
    "v * w \n",
    "= U^H g_w U  w\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4.09862852, 4.68214615, 4.70242855, 4.15529448, 3.82105494,\n",
       "       4.52767695, 4.31882443, 4.57647736, 4.95389019, 4.91595713,\n",
       "       4.35010702, 3.81069648, 4.29596588, 5.10180936, 4.27236417,\n",
       "       4.74873021, 4.94188403, 4.08582287, 5.03542755, 4.42217054])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U = matrix_U(N)\n",
    "np.abs((U.conj().transpose() / N) @ ((U @ v) * (U @ w)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4.09862852, 4.68214615, 4.70242855, 4.15529448, 3.82105494,\n",
       "       4.52767695, 4.31882443, 4.57647736, 4.95389019, 4.91595713,\n",
       "       4.35010702, 3.81069648, 4.29596588, 5.10180936, 4.27236417,\n",
       "       4.74873021, 4.94188403, 4.08582287, 5.03542755, 4.42217054])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "conv(v, w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convolution on graphs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Plan**:\n",
    "    - Define the graph Laplacian\n",
    "    - Compute the spectrum\n",
    "    - Define a Fourier transform\n",
    "    - Define convolution on a graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** From now on $G = (V, E)$ is an undirected, unweighted, simple graph."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Graph Laplacian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adjacency matrix\n",
    "\\begin{align*}\n",
    "A_{ij} = \\left\\{\n",
    "    \\begin{array}{ll}\n",
    "    1 &\\text{ if } e_{ij}\\in E\\\\\n",
    "    0 &\\text{ if } e_{ij}\\notin E\n",
    "    \\end{array}\n",
    "    \\right.\n",
    "\\end{align*}\n",
    "\n",
    "Degree matrix\n",
    "\\begin{align*}\n",
    "D_{ij} = \\left\\{\n",
    "    \\begin{array}{ll}\n",
    "    \\mbox{deg}(v_i) &\\text{ if } i=j\\\\\n",
    "    0 &\\text{ if } i\\neq j\n",
    "    \\end{array}\n",
    "    \\right.\n",
    "\\end{align*}\n",
    "\n",
    "Laplacian\n",
    "\\begin{align*}\n",
    "L &= D - A.\n",
    "\\end{align*}\n",
    "\n",
    "Normalized Laplacian\n",
    "\\begin{align*}\n",
    "L &= I - D^{-1/2} A D^{-1/2}.\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Graph spectrum, Fourier transform, and convolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Spectral decomposition of the Laplacian:\n",
    "\\begin{align*}\n",
    "L = U \\Lambda U^T\\\\\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "2. Fourier transform: if $v$ is a vector of features on the graph, then\n",
    "\\begin{align*}\n",
    "\\mathcal F (v) = U \\cdot v, \\;\\;\\mathcal F^{-1} (v) = U^T \\cdot v\\\\\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "3. Convolution with a filter $U \\cdot w$\n",
    "\\begin{align*}\n",
    "v * w = U ((U^T  w) \\odot (U^T  v) )\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "Or $g_w = \\mbox{diag}(U^T w)$ is a filter, then\n",
    "\\begin{align*}\n",
    "v * w = U g_w U^T  v\n",
    "\\end{align*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spectral-convolutional layers in PyTorch Geometric"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Problem:** Computing the spectrum is a global and very expensive property.\n",
    "\n",
    "**Goal:** Implementation as message passing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ChebConv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Original [paper](https://arxiv.org/pdf/1606.09375.pdf)\n",
    "- PyTorch [doc](https://pytorch-geometric.readthedocs.io/en/latest/modules/nn.html#torch_geometric.nn.conv.ChebConv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Goal: \n",
    "Compute $U g_w U^T x$ with $g_w = g_w(\\Lambda)$ a filter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Chebyshev approximation\n",
    "\n",
    "Chebyshev polynomials $T_k$:\n",
    "\\begin{align*}\n",
    "T_{k}(x) = 2 x T_{k-1}(x) - T_{k-2}(x), \\;\\; T_0(x) = 1, T_1(x) = x\n",
    "\\end{align*}\n",
    "\n",
    "#### Chebyshev approximation of the filter\n",
    "Aproximation of the filter:\n",
    "\\begin{align*}\n",
    "g_w(\\Lambda) = \\sum_{k=0}^K \\theta_k T_k(\\tilde \\Lambda),\\;\\;\\;\\;\\tilde \\Lambda = \\frac{2}{\\lambda_\\max} \\Lambda - I \\cdot \\lambda_\\max\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "#### Property\n",
    "If $L = U \\Lambda U^T$ then $T_k(L) = U T_k(\\Lambda) U^T$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fast approximated convolution \n",
    "\\begin{align*}\n",
    "v * w &= U g_w U^T x\n",
    "= U \\left(\\sum_{k=0}^K \\theta_k T_k(\\tilde \\Lambda) \\right)U^T x\n",
    "=\\sum_{k=0}^K  \\theta_k U  T_k(\\tilde \\Lambda) U^T x\\\\ \n",
    "&=\\sum_{k=0}^K  \\theta_k T_k(\\tilde L) x \n",
    "\\end{align*}\n",
    "\n",
    "\\begin{align*}\n",
    "\\tilde L = \\frac{2}{\\lambda_\\max} L - I\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Properties:\n",
    "- Depends on $L$ and $\\lambda_\\max$, not on $U, \\Sigma$\n",
    "- Uses only $K$-powers $\\Rightarrow$ only the $K$-th neighborhood of each node, localized filter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**As message passing:**\n",
    "![](fig/cheb_init.png)\n",
    "\n",
    "![](fig/cheb_norm.png)\n",
    "\n",
    "![](fig/cheb_forward.png)\n",
    "\n",
    "![](fig/cheb_message.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GCNConv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Original [paper](https://arxiv.org/pdf/1609.02907.pdf)\n",
    "- PyTorch [doc](https://pytorch-geometric.readthedocs.io/en/latest/modules/nn.html#torch_geometric.nn.conv.GCNConv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start from `ChebConv` and assume \n",
    "1. $K=1$ (linear approximation) so\n",
    "\\begin{align*}\n",
    "v * w \n",
    "&=\\sum_{k=0}^1  \\theta_k T_k(\\tilde L) x\n",
    "= \\theta_0 x + \\theta_1 \\tilde L x\\\\\n",
    "\\end{align*}\n",
    "\n",
    "2. $\\lambda_\\max =2$ so\n",
    "\\begin{align*}\n",
    "v * w \n",
    "&= \\theta_0 x + \\theta_1 (L - I) x\\\\\n",
    "&= \\theta_0 x - \\theta_1 D^{-1/2} A D^{1/2} x\\\\\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "3. $\\theta_0=-\\theta_1= \\theta$ so \n",
    "\\begin{align*}\n",
    "v * w = \\left(I + D^{-1/2} A D^{1/2}\\right) x \\theta\n",
    "\\end{align*}\n",
    "\n",
    "4. Renormalization of $\\theta$ by using \n",
    "\\begin{align*}\n",
    "\\tilde A&:= I + A\\\\\n",
    "\\tilde D_{ii}&:= \\sum_j \\tilde A_{ij}\n",
    "\\end{align*}\n",
    "so \n",
    "\\begin{align*}\n",
    "v * w = \\left(D^{-1/2} A D^{1/2}\\right) x \\theta\n",
    "\\end{align*}\n",
    "\n",
    "If $x$ is a $F$-dimensional feature vector, and we want an $F'$-dimensional feature vector as output:\n",
    "use $W'\\in \\mathbb{R}^{F\\times F'}$\n",
    "\\begin{align*}\n",
    "v * w = \\left(D^{-1/2} A D^{1/2}\\right) x \\Theta\n",
    "\\end{align*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nodewise:\n",
    "    ![image.png](fig/gcn_nodewise.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### As message passing\n",
    "See Tutorial3"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
